This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
The objective of this notebook is to discretize HRUs in a watershed and create a weighting matrix that describes the hillslope-stream network connectivity in a catchment. This notebook is meant to replace the spatial discretization workflow (DynatopSpatialFunctionExplicitReaches) for the modified dynamic TOPMODEL described by Mahoney et al. (2022) J. Hydrol.
To run this notebook, it is necessary to first install kernels for Python and R and install arcpy to the environment from which the notebook is run (see instructions here). Instances when switching the kernel is necessary will be denoted in a Markdown cell.
We will be using both R and python in this example. Python will be called via the reticulate package. We will be using the python installation from arcgis pro. A major advantage of using the R Notebook document over, for example, Jupyter Notebook is the ability to store variables generated by both python and R in the environment. Thus we can use variables generated in R with the arcgis API, have the user-friendliness of Rstudio, and have the replicability of a Notebook document.
setwd('C:/Users/david/OneDrive/Desktop/Mahoney Research/2-Mahoney-EPA/modeling-streamflow-permanence/')
# Load in libraries
library(reticulate)
library(sp)
library(raster)
library(dynatopmodel)
library(topmodel)
library(ecbtools)
library(cluster)
library(factoextra)
library(tidyverse)
library(NbClust)
library(basicClEval)
library(ggplot2)
library(mapview)
library(doParallel)
library(foreach)
use_python('C:/Program Files/ArcGIS/Pro/bin/Python/envs/arcgispro-py3/python.exe',required=T)
Warning in Sys.setlocale("LC_CTYPE", ctype) :
using locale code page other than 65001 ("UTF-8") may cause problems
arcpy <- import("arcpy")
arcgis <- import("arcgis")
py_math <- import("math")
py_numpy <- import("numpy")
ipython <- import("IPython")
Note: This code will calculate the TWI using the D8 flow direction. It is possible to import a DINF or MFD raster here as well and use this to run the analysis, assuming that it is already in the GDB.
# Read in the rasters that we'd need for HRU creation
dem <- raster('SpatialInputData/fr1meterDEM.tif')
soils <- raster('SpatialInputData/FRSoils1mMask.tif')
TWI.dinf <- raster('C:/Users/david/OneDrive/Desktop/Mahoney Research/2-Mahoney-EPA/modeling-streamflow-permanence/SpatialInputData/fr_twi_dinf.tif')
TWI.dinf[TWI.dinf==0] <- NA # We do this because the areas outside of the raster have a value 0 and we want these to be set to NA
Now, we will switch over to R and create an HRU map of for the watershed. Make sure to switch the kernel from python (esri) to R
Relevant libraries include spatial and raster libraries, cluster analyses, and tidyverse. Note, we are creating HRUs with the following rasters:
DEM
Soils
TWI
Other relevant layers that could be used include:
LULC
Ksat
Porosity
Aspect
Curvature
Temperature
Precipitation
Note: there is some degree of customizability here. Any relevant layer could be added.
The total within sum of squares for each cluster is used to estimate the variability for each iteration.
Note: this code takes a decent amount of time to run.
wss.raster <- function(k) {
# Run the K-means analysis
kmeans.raster <- raster.kmeans(raster.stack, k, iter.max=10, nstart = 10, geo = T, geo.weight = 1)
# Load in the data and get into vector form
HRU.vector <- data.frame('hru' = as.vector(kmeans.raster))
TWI.vector <- data.frame('TWI' = as.vector(TWI.dinf.cut))
soils.vector <- data.frame('soils' = as.vector(soils))
elev.bands.vector <- data.frame('Elev.bands' = as.vector(elev.bands))
# Put the raster vectors into a cleaned data frame
rasters.df <- data.frame(HRU.vector,TWI.vector,soils.vector,elev.bands.vector)
rasters.df.clean <- rasters.df[complete.cases(rasters.df),]
# Calculate the within cluster sum of squares
within.ss <- wcss(rasters.df.clean[,-1],rasters.df.clean[,1])
# Return the total within cluster sum of squares
return(within.ss$WCSS)
}
set.seed(123)
# Compute and plot within sum squares for k = 1 to k = 40
k.values <- 1:30
# extract within sum of squares for 1-40 clusters
wss_values <- map_dbl(k.values,wss.raster)
Warning: Quick-TRANSfer stage steps exceeded maximum (= 20786250)
Warning: Quick-TRANSfer stage steps exceeded maximum (= 20786250)
# plot within sum of squares vs cluster (HRU) number
plot(k.values,wss_values,type='b',pch=19,frame=F,xlab='Number of clusters',ylab='Total within-cluster sum of sq')
We will look at the various clusters created and how they appear within relation to one another
# Load in the data and get into vector form
HRU.vector <- data.frame('hru' = as.vector(HRU))
TWI.vector <- data.frame('TWI' = as.vector(TWI.dinf.cut))
TWI.nc.vector <- data.frame('TWI.nc'=as.vector(TWI.dinf))
soils.vector <- data.frame('soils' = as.vector(soils))
elev.bands.vector <- data.frame('Elev.bands' = as.vector(elev.bands))
elev.bands.nc.vector <- data.frame('Elev.bands.nc' = as.vector(dem))
# Put the raster vectors into a cleaned data frame
rasters.df <- data.frame(HRU.vector,TWI.vector,TWI.nc.vector,soils.vector,elev.bands.vector,elev.bands.nc.vector)
rasters.df.clean <- rasters.df[complete.cases(rasters.df),]
# Plot the elev bands vs twi with HRU as color
ggplot(data=rasters.df.clean, aes(x=TWI.nc,y=Elev.bands.nc,col=hru))+geom_point()
# Note - the straight lines are due to the cuts that we creted. We could just as easily not include these cuts and other geometries would come through
In this step, we will first sort HRUs by TWI value (high TWI = high number HRU), then we will “burn in” a stream network that has been previously delineated. This represents the geomorphic stream channel, to which the watershed uplands will contribute water.
The raster need to have the following characteristics:
The resolution, extent, and snapping should be the same as the DEM
The raster is derived from a Flow accumulation raster with an upstream contributing area threshold applied
The threshold should create a stream network that is similar to the stream mapped in the field
The reaches should be created using the stream link/order tools in ArcGIS
The shapefile created from these stream tools should be buffered such that a contiguous stream network raster can be created
Here, a buffer of 2 cell widths was used (~1.524 * 2 m)
## Calculate the mean TWI for each HRU - relabel the HRUs such that high valued TWI is the highest number HRU and low valued TWI is the lowest HRU
## Note - this isn't necessary for dynatopmodel, but may be necessary for dynatop
mean.twi <- data.frame(matrix(nrow=length(unique(HRU)),ncol=5))
names(mean.twi) <- c('HRU.name','atb.bar.unsort','id.new','atb.bar','id')
mean.twi$HRU.name <- unique(HRU)
hillslope.hru.matrix <- as.matrix(HRU)
TWI.dinf.matrix <- as.matrix(TWI.dinf)
mean.twi$atb.bar.unsort <- zonal(x=TWI.dinf,z=HRU,fun='mean')[,2]
mean.twi$atb.bar <- sort(mean.twi$atb.bar.unsort,decreasing=F)
mean.twi$id.new <- match(mean.twi$atb.bar.unsort,mean.twi$atb.bar)
mean.twi$id <- 1:length(unique(HRU))
reclassify.matrix <- data.frame(matrix(nrow=length(unique(HRU)),ncol=2))
names(reclassify.matrix) <- c('is','becomes')
reclassify.matrix$is <- mean.twi$HRU.name
reclassify.matrix$becomes <- mean.twi$id.new
HRU <- reclassify(HRU,reclassify.matrix)
## Overlay the stream network atop the Hru raster
stream <- raster('SpatialInputData/FR_stream_network.tif')
stream[stream==0] <- NA
HRU <- HRU*10+max(max(unique(stream)),100) # Uniquely identifying the upland HRUs
HRU[stream>0] <- stream[stream>0] # Identify the cells from the HRU raster that overlap with the stream, set those equal to the reach number
writeRaster(HRU,"SpatialInputData/FR_HRU.tif",overwrite=T)
image(HRU)
twi.zonal <- zonal(x=TWI.dinf,z=HRU,fun='mean')
TWI.zonal <- reclassify(HRU,twi.zonal)
image(TWI.zonal)
stream.lump <- stream
stream.lump[stream >= 1] <- 1
HRU.lump <- HRU
HRU.lump[stream.lump>0] <- stream.lump[stream.lump>0]
writeRaster(HRU.lump,"SpatialInputData/FR_HRU_lump.tif",overwrite=T)
writeRaster(stream.lump,"SpatialInputData/FR_stream_lump.tif",overwrite=T)
image(HRU.lump)
Note - I’m going to do this with the mapview function in R, but this could be done using the arcpy and arcgis modules in python. The reason for not using python is because if we interrupt the kernel we’ll have to reload variables from above
mapView(HRU, maxpixels = 2102000) + mapView(TWI.zonal, maxpixels = 2102000) +
mapView(HRU.lump,maxpixels = 2102000) + mapView(stream.lump, maxpixels = 2102000) +
mapView(stream,maxpixels = 2102000)
Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs
Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
Discarded datum World Geodetic System 1984 in Proj4 definition
In this step, we will create the weighting matrix to run the dynatopmodel analysis. We do this twice. Once for the “Explicit matrix” which shows the connectivity of upland cells that contribute to a specific reach. A second time for the “Lumped Matrix” which shows the connectivity of cells to a “lumped” stream network. As of 8/2/2022, the “dynatopmodel” package is not compatible with explicit reaches.
This process returns the value of the HRU immediately downstream of a cell. This algorithm strictly works for a D8 flow direction, but could be expanded for DInf at a later time.
Make sure to switch over to the R kernel
# First, create the explicit matrix
library(sp)
library(raster)
HRU <- raster('SpatialInputData/FR_HRU.tif')
fdr <- raster('SpatialInputData/fr_fdr')
HRU.matrix <- as.matrix(HRU)
fdr.matrix <- as.matrix(fdr)
fdr.downstream.matrix <- matrix(nrow=nrow(HRU.matrix),ncol=ncol(HRU.matrix))
for (i in 1:nrow(HRU.matrix)) {
for (j in 1:ncol(HRU.matrix)) {
if (is.na(fdr.matrix[i,j])) {
fdr.downstream.matrix[i,j] <- NA
} else {
if (fdr.matrix[i,j]==1) {
fdr.downstream.matrix[i,j] <- HRU.matrix[i,j+1]
} else if (fdr.matrix[i,j]==2) {
fdr.downstream.matrix[i,j] <- HRU.matrix[i+1,j+1]
} else if (fdr.matrix[i,j]==4) {
fdr.downstream.matrix[i,j] <- HRU.matrix[i+1,j]
} else if (fdr.matrix[i,j]==8) {
fdr.downstream.matrix[i,j] <- HRU.matrix[i+1,j-1]
} else if (fdr.matrix[i,j]==16) {
fdr.downstream.matrix[i,j] <- HRU.matrix[i,j-1]
} else if (fdr.matrix[i,j]==32) {
fdr.downstream.matrix[i,j] <- HRU.matrix[i-1,j-1]
} else if (fdr.matrix[i,j]==64) {
fdr.downstream.matrix[i,j] <- HRU.matrix[i-1,j]
} else if (fdr.matrix[i,j]==128) {
fdr.downstream.matrix[i,j] <- HRU.matrix[i-1,j+1]
} else {
fdr.downstream.matrix[i,j] <- NA
}
}
}
}
In this step, we will create a weighting matrix that outputs the proportion of flow generated in one cell that will be redistributed to downstream cells.
This works by:
identifying all cells that belong to a certain HRU
identifying all cells which belong to a downstream HRU
identifying the cells where both of these conditions are true
returning the number of cells where both are true
In this regard, the output will be all cells which belong to HRU i that flow into HRU j.
Note: I rewrote some of this code. The run time goes from about 50 minutes to approximately 50 seconds.
print('Running weighting matrix...')
[1] "Running weighting matrix..."
no.HRU.explicit <- length(unique(HRU))
HRU.explicit.list <- unique(HRU)
downstream.weighting.matrix <- matrix(nrow=no.HRU.explicit,ncol=no.HRU.explicit)
start.time <- Sys.time()
for (i in 1:no.HRU.explicit) {
print(HRU.explicit.list[i])
for (j in 1:no.HRU.explicit) {
downstream.weighting.matrix[i,j] <- sum((HRU.matrix==HRU.explicit.list[i])&(fdr.downstream.matrix==HRU.explicit.list[j]),na.rm=T)
}
}
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14
[1] 15
[1] 16
[1] 17
[1] 18
[1] 19
[1] 20
[1] 21
[1] 22
[1] 23
[1] 24
[1] 25
[1] 26
[1] 27
[1] 28
[1] 29
[1] 30
[1] 31
[1] 32
[1] 33
[1] 110
[1] 120
[1] 130
[1] 140
[1] 150
[1] 160
[1] 170
[1] 180
[1] 190
[1] 200
[1] 210
[1] 220
[1] 230
[1] 240
[1] 250
[1] 260
[1] 270
[1] 280
[1] 290
[1] 300
end.time <- Sys.time()
end.time-start.time
Time difference of 53.847 secs
options("scipen"=100, "digits"=5)
rownames(downstream.weighting.matrix) <- unique(HRU)
colnames(downstream.weighting.matrix) <- unique(HRU)
weighting.matrix <- downstream.weighting.matrix/rowSums(downstream.weighting.matrix)
rownames(weighting.matrix) <- unique(HRU)
colnames(weighting.matrix) <- unique(HRU)
rowSums(weighting.matrix)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
21 22 23 24 25 26 27 28 29 30 31 32 33 110 120 130 140 150 160 170
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
180 190 200 210 220 230 240 250 260 270 280 290 300
1 1 1 1 1 1 1 1 1 1 1 1 1
Note that the current implementation of dynatopmodel does not support the explicit streams, as far as I am aware. Dynatop may, but I haven’t investigated the model in depth as of 8/3/2022
# Next, create the lumped matrix
HRU.lump <- raster('SpatialInputData/FR_HRU_lump.tif')
fdr <- raster('SpatialInputData/fr_fdr')
HRU.lump.matrix <- as.matrix(HRU.lump)
fdr.matrix <- as.matrix(fdr)
fdr.downstream.lump.matrix <- matrix(nrow=nrow(HRU.lump.matrix),ncol=ncol(HRU.lump.matrix))
for (i in 1:nrow(HRU.lump.matrix)) {
for (j in 1:ncol(HRU.lump.matrix)) {
if (is.na(fdr.matrix[i,j])) {
fdr.downstream.lump.matrix[i,j] <- NA
} else {
if (fdr.matrix[i,j]==1) {
fdr.downstream.lump.matrix[i,j] <- HRU.lump.matrix[i,j+1]
} else if (fdr.matrix[i,j]==2) {
fdr.downstream.lump.matrix[i,j] <- HRU.lump.matrix[i+1,j+1]
} else if (fdr.matrix[i,j]==4) {
fdr.downstream.lump.matrix[i,j] <- HRU.lump.matrix[i+1,j]
} else if (fdr.matrix[i,j]==8) {
fdr.downstream.lump.matrix[i,j] <- HRU.lump.matrix[i+1,j-1]
} else if (fdr.matrix[i,j]==16) {
fdr.downstream.lump.matrix[i,j] <- HRU.lump.matrix[i,j-1]
} else if (fdr.matrix[i,j]==32) {
fdr.downstream.lump.matrix[i,j] <- HRU.lump.matrix[i-1,j-1]
} else if (fdr.matrix[i,j]==64) {
fdr.downstream.lump.matrix[i,j] <- HRU.lump.matrix[i-1,j]
} else if (fdr.matrix[i,j]==128) {
fdr.downstream.lump.matrix[i,j] <- HRU.lump.matrix[i-1,j+1]
} else {
fdr.downstream.lump.matrix[i,j] <- NA
}
}
}
}
Note: I have rewritten the code here. The previous run time for this chunk was 10 minutes. Now it is about 10 seconds.
print('Running weighting matrix...')
[1] "Running weighting matrix..."
no.HRU.lump <- length(unique(HRU.lump))
list.HRU.lump <- unique(HRU.lump)
downstream.weighting.matrix.lump <- matrix(nrow=length(unique(HRU.lump)),ncol=length(unique(HRU.lump)))
start.time <- Sys.time()
for (i in 1:no.HRU.lump) {
print(list.HRU.lump[i])
for (j in 1:no.HRU.lump) {
downstream.weighting.matrix.lump[i,j] <- sum((HRU.lump.matrix==list.HRU.lump[i])&(fdr.downstream.lump.matrix==list.HRU.lump[j]),
na.rm=T)
}
}
[1] 1
[1] 110
[1] 120
[1] 130
[1] 140
[1] 150
[1] 160
[1] 170
[1] 180
[1] 190
[1] 200
[1] 210
[1] 220
[1] 230
[1] 240
[1] 250
[1] 260
[1] 270
[1] 280
[1] 290
[1] 300
end.time <- Sys.time()
end.time-start.time
Time difference of 10.993 secs
options("scipen"=100, "digits"=5)
rownames(downstream.weighting.matrix.lump) <- unique(HRU.lump)
colnames(downstream.weighting.matrix.lump) <- unique(HRU.lump)
weighting.matrix.lump <- downstream.weighting.matrix.lump/rowSums(downstream.weighting.matrix.lump)
rownames(weighting.matrix.lump) <- unique(HRU.lump)
colnames(weighting.matrix.lump) <- unique(HRU.lump)
rowSums(weighting.matrix.lump)
1 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
300
1
Next, we will create a “Groups” matrix which will specify the initial parameters for dynamic TOPMODEL.
Note - I did notice some differences between atb.bar from the 20 HRUs discretized here and the 2- HRUs discretized with dynatopmodel. It is related to how the HRUs are discretized. There is more variability with the dynatopmodel approach.
Items we will need to include in the table:
id - the ID of the channel, typically channels are ranked 0-100 and HRUs are ranked based on TWI with ID > 100
tag - same as id
chan.no - the channel number (note - only for IDs classified as channels)
area_pc - the percent area that the HRU takes up
area - the area m^2 that the HRU takes up
sbar - the average slope of the HRU (m/m)
atb.bar - the average TWI for the HRU
gauge.id - the rain gauge contributing to the HRU
catch.id - the catchment in which the HRU is situated
srz_max - maximum depth of root zone parameter value initialized as 0.1
ln_t0 - ln t0 parameter value initialized as 7
m - exponential decline of transmissivity parameter initialized as 0.01
srz0 - initial root zone storage parameter initialized as 0
td - td parameter initialized as 1
vchan - v chan parameter (for channels only) parameterized as 1000
vof - overland v parameter (for HRUs only) parameterized as 100
k0 - k0 parameter initialized as 1e+8
CD - CD parameter (unused) initialized as 0.1
sd_max - sd_max parameter initialized as 0.5
pe_fact - potential evapotranspiration factor used to scale PET, initialized as 1
vof_fact - vof factor initialized as 1
rain_fact - rain factor used to scale precip amount, initialized as 1
mann.n - mannings n initialized as 0.01
S0 - initial slope, parameterized as 0.1
ex_max - ex max parameterized as 1
# Enter the cell size
cell.size <- as.numeric(readline('Enter the resolution of the raster in m (for 10-m raster enter 10): '))
1.524
# Create functions to determine the area of the raster
area_cell <- function(raster.matrix,raster_ID) {
area_c <- sum(raster.matrix==raster_ID, na.rm=T)
return(area_c)
}
# Create function to determine the zonal mean (unused)
mean_raster <- function(raster.matrix, HRU.matrix, raster_ID) {
mean_rast <- mean(raster.matrix[HRU.matrix==raster_ID],na.rm=T)
return(mean_rast)
}
# Calculate the slope
slope <- terrain(dem,opt='slope',units='tangent')
# Initial ize the groups
groups.explicit <- data.frame(matrix(nrow=length(unique(HRU)),ncol=26))
names(groups.explicit) <- c('id','tag','chan.no','order','area_pc','area','sbar','atb.bar','gauge.id','catch.id','srz_max','ln_t0','m','srz0','td','vchan','vof','k0','CD','sd_max','pe_fact','vof_fact','rain_fact','mann.n','S0','ex_max')
# Assign the groups
groups.explicit$id <- unique(HRU)
groups.explicit$tag <- unique(HRU)
groups.explicit$chan.no[1:length(unique(stream))] <- unique(stream)
groups.explicit$order <- 1:length(unique(HRU))
groups.explicit$area_pc <- sapply(X=groups.explicit$id,FUN=area_cell,raster.matrix=HRU.matrix)/sum(!is.na(HRU.matrix))*100
groups.explicit$area <- groups.explicit$area_pc*sum(!is.na(HRU.matrix))*cell.size^2/100
groups.explicit$sbar <- zonal(x=slope,z=HRU,fun='mean')[,2]
groups.explicit$atb.bar[(1+length(unique(stream))):length(unique(HRU))] <- zonal(x=TWI.dinf,z=HRU,fun='mean')[(1+length(unique(stream))):length(unique(HRU)),2]
groups.explicit$gauge.id <- 1
groups.explicit$catch.id <- 1
groups.explicit$srz_max <- 0.1
groups.explicit$ln_t0 <- 7
groups.explicit$m <- 0.01
groups.explicit$srz0 <- 0
groups.explicit$td <- 1
groups.explicit$vchan[1:length(unique(stream))] <- 1000
groups.explicit$vof[(1+length(unique(stream))):length(unique(HRU))] <- 100
groups.explicit$k0 <- 1e+8
groups.explicit$CD <- 0.1
groups.explicit$sd_max <- 0.5
groups.explicit$pe_fact <- 1
groups.explicit$vof_fact <- 1
groups.explicit$rain_fact <- 1
groups.explicit$mann.n <- 0.01
groups.explicit$S0 <- 0.1
groups.explicit$ex_max <- 1
# Calculate the slope
slope <- terrain(dem,opt='slope',units='tangent')
# Initialize the groups
groups.lump <- data.frame(matrix(nrow=no.HRU.lump,ncol=26))
names(groups.lump) <- c('id','tag','chan.no','order','area_pc','area','sbar','atb.bar','gauge.id','catch.id','srz_max','ln_t0','m','srz0','td','vchan','vof','k0','CD','sd_max','pe_fact','vof_fact','rain_fact','mann.n','S0','ex_max')
# Assign the groups
groups.lump$id <- list.HRU.lump
groups.lump$tag <- list.HRU.lump
groups.lump$chan.no[1:length(unique(stream.lump))] <- unique(stream.lump)
groups.lump$order <- 1:no.HRU.lump
groups.lump$area_pc <- sapply(X=groups.lump$id,
FUN=area_cell,raster.matrix=HRU.lump.matrix)/sum(!is.na(HRU.lump.matrix))*100
groups.lump$area <- groups.lump$area_pc*sum(!is.na(HRU.lump.matrix))*cell.size^2/100
groups.lump$sbar <- zonal(x=slope,z=HRU.lump,fun='mean')[,2]
groups.lump$atb.bar[(1+length(unique(stream.lump))):length(unique(HRU.lump))] <- zonal(x=TWI.dinf,z=HRU.lump,fun='mean')[(1+length(unique(stream.lump))):length(unique(HRU.lump)),2]
groups.lump$gauge.id <- 1
groups.lump$catch.id <- 1
groups.lump$srz_max <- 0.1
groups.lump$ln_t0 <- 7
groups.lump$m <- 0.01
groups.lump$srz0 <- 0
groups.lump$td <- 1
groups.lump$vchan[1:length(unique(stream.lump))] <- 1000
groups.lump$vof[(1+length(unique(stream.lump))):length(unique(HRU.lump))] <- 100
groups.lump$k0 <- 1e+8
groups.lump$CD <- 0.1
groups.lump$sd_max <- 0.5
groups.lump$pe_fact <- 1
groups.lump$vof_fact <- 1
groups.lump$rain_fact <- 1
groups.lump$mann.n <- 0.01
groups.lump$S0 <- 0.1
groups.lump$ex_max <- 1
In this step, we will create routing matrix that bins lists the length of each cell to the outlet and the frequency.
# Determine the flow length from each cell in the stream to the catchment outlet
outFlowLength <- arcpy$sa$FlowLength('SpatialInputData/fr_fdr',"DOWNSTREAM")
# Note - the next line only needs to be saved once.
#outFlowLength$save(paste0(getwd(),'/SpatialInputData/Flow_length.tif'))
Flow.length <- raster('SpatialInputData/Flow_length.tif')
Flow.length.stream <- stream.lump
Flow.length.stream[!is.na(stream.lump)] <- Flow.length[stream.lump==1]
writeRaster(Flow.length.stream,'SpatialInputData/Flow_length_stream.tif', overwrite =T)
input.units <- readline('Was the original DEM in units of ft or m (enter ft or m): ')
ft
breaks.hist <- 10 # User defined number of breaks
RoutingTable <- data.frame(matrix(nrow=breaks.hist,ncol=2))
names(RoutingTable) <- c('flow.len', 'prop')
if(input.units == 'ft') {
Flow.length.stream.fix <- Flow.length.stream * 0.3048
Flow.length.vector <- na.omit(as.numeric(as.vector(Flow.length.stream.fix)))
RoutingHist <- hist(Flow.length.vector,breaks=seq(min(Flow.length.vector),max(Flow.length.vector),length.out=breaks.hist+1))
RoutingTable$flow.len <- RoutingHist$breaks[-1]
RoutingTable$prop <- RoutingHist$counts/sum(RoutingHist$counts)
} else if (input.units == 'm') {
Flow.length.vector <- na.omit(as.numeric(as.vector(Flow.length.stream)))
RoutingHist <- hist(Flow.length.vector,breaks=seq(min(Flow.length.vector),max(Flow.length.vector),length.out=breaks.hist+1))
RoutingTable$flow.len <- RoutingHist$breaks[-1]
RoutingTable$prop <- RoutingHist$counts/sum(RoutingHist$counts)
} else {
'invalid response, RoutingTable not calculated.'
}
We will create a list that compThe original DynatopSpatialFunctionExplicitReaches function outputs the following information:
DEM
Soils
layers - contains all layers to make HRUs
DEM
Soils
TWI (named atb)
disc - contains all elements to make the weighting matrix and groups table from the lumped parameterization
groups - the groups matrix
layers - same as above
chans - multiband raster with 1) channel labels (just equal to 1 for lumped) and 2) channel props, which describes the proportions of the cell which is occupied by the channel (here this can just be set equal to one if a 1.524-m raster is being used).
cuts - the number of cuts made to the layers
area.thresh - the area threshold below which HRUs will be lumped with other HRUs. This isn’t utilized here and can be set equal to 0.
hru - the HRU raster generated from the discretize process
weights - the weighting matrix
RoutingTable - a routing table showing the flow distance to the outlet and the proportion of channel cells that fall within the distance bin - this had 30 cuts
explicit.disc - contains all elements to make the weighting matrix and groups table from the explicit parameterization
groups - explicit groups matrix
layers - should be more or less the same as disc
chans - This time should be explicit with respect to Reach ID
cuts - should be the same as disc
area.thresh - should be same as disc
hru - This will be different than disc bc there should be explicit representations for each reach
weights - explicit weighting matrix
explicit.RoutingTable - same as above, this had 5 cuts, however
explicit.ChanTable - this is a table containing all of the reach data, including: “link_no,” “ds_link,” “us_link_1,” “us_link_2,” “stream_order,” “stream_length_ft,” “stream_drop_ft,” “stream_slope_ftpft,” “stream_straightline_length_ft,” “US_area_m2,” “width_m2”
# DEM
DEM <- dem
# Soils
Soils <- soils
# Layers
layers <- raster.stack
# disc
disc <- list()
disc$groups <- groups.lump
disc$layers <- layers
disc$chans <- addLayer(stream.lump,stream.lump)
names(disc$chans) <- c('chans','chanprops')
disc$cuts <- data.frame(matrix(nrow=1,ncol=length(names(layers))))
colnames(disc$cuts) <- names(layers)
for (i in names(disc$cuts)) {print(i);disc$cuts[1,i] <- length(unique(layers[[i]]))}
[1] "TWI.dinf.cut"
[1] "soils"
[1] "elev.bands"
disc$cuts <- as.matrix(disc$cuts)
disc$area.thresh <- 0
disc$hru <- HRU.lump
weighting.matrix.lump.in <- weighting.matrix.lump
colnames(weighting.matrix.lump.in)[1] <- 'R'
rownames(weighting.matrix.lump.in)[1] <- 'R'
disc$weights <- as.matrix(weighting.matrix.lump.in)
# RoutingTable
RoutingTable <- RoutingTable
# explicit.disc
explicit.disc <- list()
explicit.disc$groups <- groups.explicit
explicit.disc$layers <- layers
explicit.disc$chans <- addLayer(stream,stream.lump)
names(explicit.disc$chans) <- c('chans','chanprops')
explicit.disc$cuts <- disc$cuts
explicit.disc$area.thresh <- 0
explicit.disc$hru <- HRU
weighting.matrix.explicit.in <- weighting.matrix
names.explicit.matrix <- paste0('R',unique(stream))
rownames(weighting.matrix.explicit.in)[1:length(unique(stream))] <- names.explicit.matrix
colnames(weighting.matrix.explicit.in)[1:length(unique(stream))] <- names.explicit.matrix
explicit.disc$weights <- as.matrix(weighting.matrix.explicit.in)
# explicit.RoutingTable (unused)
explicit.RoutingTable <- RoutingTable
# explicit.ChanTable
drn <- shapefile('SpatialInputData/FR4000StreamNet.shp')
bkf.width <- 12.16*(drn$USAreaM2/1000^2*0.386102)^0.42*0.3048 # estimate the width of the channel using the regional curves from Ashland Berry's thesis. Note this isn't perfect bc this equation is being extrapolated to smaller channels.
explicit.ChanTable <- as.data.frame(cbind('link_no'= drn$LINKNO,'ds_link'= drn$DSLINKNO,
'us_link_1' = drn$USLINKNO1,'us_link_2'=drn$USLINKNO2,
'stream_order'= drn$strmOrder,
'stream_length_m' =drn$length_m,
'stream_length_ft'=drn$Length,
'stream_drop_ft'=drn$strmDrop,
'stream_slope_ftpft'=drn$Slope,
'stream_straightline_length_ft'=drn$StraightL,
'US_area_m2'=drn$USAreaM2,
'width_m2' = bkf.width))
# Compile the final spatial table
dynatop.spatial <- list('DEM'=DEM,'Soils'=Soils,'layers'=layers,'disc'=disc,
'RoutingTable'=RoutingTable,'explicit.disc'=explicit.disc,
'explicit.RoutingTable'=explicit.RoutingTable,
'explicit.ChanTable'=explicit.ChanTable)
This chunk is a scratch workspace for old code
tot.cells <- data.frame(matrix(nrow=no.HRU.lump,ncol=2))
names(tot.cells) <- c('tot.cells','tot.cells.2')
for (j in 1:no.HRU.lump) {
print(list.HRU[j])
tot.cells$tot.cells[j] <- sum((HRU.lump.matrix==list.HRU[1])&(fdr.downstream.lump.matrix==list.HRU[j]),na.rm=T)
#tot.cells$tot.cells.2[j] <- length(which((HRU.lump.matrix==unique(HRU.lump)[1])&(fdr.downstream.matrix==unique(HRU.lump)[j])))
}
parallel::detectCores()
n.cores <- parallel::detectCores() - 1
#create the cluster
my.cluster <- parallel::makeCluster(
n.cores,
type = "PSOCK"
)
#check cluster definition (optional)
print(my.cluster)
#register it to be used by %dopar%
doParallel::registerDoParallel(cl = my.cluster)
#check if it is registered (optional)
foreach::getDoParRegistered()
get.ds.cells <- function(i,HRU.matrix,HRU,fdr.downstream.matrix) {
col <- matrix(nrow=1,length(unique(HRU)))
for (j in 1:length(unique(HRU))) {
temp.col <- length(which((HRU.matrix==unique(HRU)[i])&(fdr.downstream.matrix==unique(HRU)[j])))
col[j] <- temp.col
}
return(col)
}
start.time <- Sys.time()
downstream.weighting.matrix.test <- foreach(i=1:length(unique(HRU)),
.combine =rbind) %dopar% {
temp.col <- get.ds.cells(i,HRU.matrix,HRU,fdr.downstream.matrix)
}
end.time <- Sys.time()
end.time-start.time
stopCluster(cl = my.cluster)
options("scipen"=100, "digits"=5)
rownames(downstream.weighting.matrix.test) <- unique(HRU)
colnames(downstream.weighting.matrix.test) <- unique(HRU)
weighting.matrix.test <- downstream.weighting.matrix.test/rowSums(downstream.weighting.matrix.test)
rownames(weighting.matrix.test) <- unique(HRU)
colnames(weighting.matrix.test) <- unique(HRU)